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We propose a mechanism for binding of diatomic ligands to heme based on a dynamical orbital 
selection process. This scenario may be described as bonding determined by local valence fluctua- 
tions. We support this model using linear-scaling first-principles calculations, in combination with 
dynamical mean- field theory, applied to heme, the kernel of the hemoglobin metalloprotein central 
to human respiration. We find that variations in Hund's exchange coupling induce a reduction of 
the iron 3d density, with a concomitant increase of valence fluctuations. We discuss the comparison 
between our computed optical absorption spectra and experimental data, our picture accounting for 
the observation of optical transitions in the infrared regime, and how the Hund's coupling reduces, 
by a factor of five, the strong imbalance in the binding energies of heme with CO and O2 ligands. 



Metalloporphyrin systems, such as heme, play a cen- 
tral role in biochemistry. The ability of such molecules 
to reversibly bind small ligands is of great interest, par- 
ticularly in the case of heme which binds diatomic lig- 
ands such as oxygen and carbon monoxide. Heme acts 
as a transport molecule for oxygen in human respiration, 
while carbon monoxide inhibits this function. Despite 
intensive studies [1-3], the binding of the iron atom at 
centre of the heme molecule to O2 and CO ligands re- 
mains poorly understood. In particular, one problem 
obtained with density functional theory [4] (DFT) ap- 
proaches to ligand binding of heme is that the difference 
in the binding energy (AAE) of carboxy-heme and oxy- 
heme is very large, and the theory predicts an unrealistic 
binding affinity to CO, several orders of magnitude larger 
than to O2 [5, 6]. 

Recent progress has been made to cure this problem 
using DFT+/7 for molecular systems [7, 8], with which it 
was found that the inclusion of many body effects in the 
calculations reduced the imbalance between O2 and CO 
affinities [9]. Inclusion of conformal modifications, such 
as the Fe-C-0 binding angle [10], or the deviation of the 
Fe atom from the porphyrin plane, were also shown to 
affect CO and O2 binding energies. 

A general problem encountered by DFT is the strong 
dependence of the energetics and the spin state on small 
changes in the geometry. In particular, traditional DFT 
fails to describe the correct high-spin ground state of 
heme molecules. DFT+U provides an improved descrip- 
tion [7, 11], but is known to overestimate magnetic mo- 
ments and gives often artificial and non physical spin- 
symmetry-broken states. Moreover, the rotationally- 
invariant DFT+U methodology does not capture well the 
effect of the Hund's coupling J, which is known to be 
large in iron based systems. It was recently shown that 
the effect of strong correlations are not always driven by 



the Coulomb repulsion U alone, but in some cases act in 
combination with the Hund's coupling J [12-14]. Under- 
standing the effect of strong correlations in heme, and 
in particular how the symmetry of the highest occupied 
molecular orbital (HOMO) is affected by U and J, is 
important in the context of describing the CO binding, 
which was shown to be strongly dependent on the HOMO 
symmetry [15]. 

Recent progress has been made in this direction by dy- 
namical mean field theory [16] (DMFT), combined with 
DFT (DFT+DMFT) which can refine the description of 
the charge and spin of correlated ions, and describes in a 
remarkable way the strong correlations, induced by both 
U and J. Also, DFT can only describe a static magnetic 
moment associated with a spin symmetry broken state, 
and requires the inclusion of the spin-orbit interaction to 
explain a change of spin states [17]. This is not necessary 
at the DMFT level, which describes both static and fluc- 
tuating magnetic moments within the same framework. 

In this work, we extend the DFT+U analysis by 
means of the combination of state-of-the-art linear scal- 
ing DFT [18] with DMFT, and apply this methodology 
to heme. The methodology builds upon our earlier works 
[19] and is described in detail in the supplementary ma- 
terial. 

Although DFT+DMFT has been widely used to 
study solids, in this study we apply our real-space 
DFT+DMFT implementation to a moderately large 
molecule, extending the scope of applicability of DMFT 
to biology in an unprecedented manner. DMFT allows 
the quantum and thermal fluctuations, missing in zero- 
temperature DFT calculations, to be recovered. More- 
over, it includes within the calculation both the Coulomb 
repulsion U and the Hund's coupling J. Which of U or 
J drives the many body effects in heme [14] remains an 
open question, paramount to understanding ligand bind- 
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FIG. 1: Orbital selection scenario: Dependence of a) 
the iron 3d subspace occupancy rid and b) the effective spin 
quantum number s on the Hund's coupling J, for both unli- 
gated and ligated heme models. The physically relevant re- 
gion 0.5 eV < J < 1 eV is highlighted in yellow. Isosur faces of 
the real-space representation of the electronic spectral density 
of the HOMO of FeP-d for c) J = eV and d) J = 0.8 eV. 
The large central sphere shows the location of the iron atom, 
and the four blue spheres indicate nitrogen atoms. 

ing, that we address in this work. Methods are available 
to obtain U and J parameters appropriate to DMFT [20], 
but in this work we focus on the dependence of the results 
with the Hund's coupling J, and we verify that our calcu- 
lations are not sensitive to the Coulomb repulsion U or to 
the temperature T [21] . The key question that we address 
in this work is: to what extent does the Hund's coupling, 
so far neglected in all studies applied to heme, affect the 
binding of heme to O2 and CO ligands, and in particular 
does J reduce the strong affinity for CO binding? If not 
specified otherwise, we use a similar value U = 4 eV to 
those previously computed for DFT+U [7], and ambient 
temperature T = 294K. The methodology is described 
in detail the supplementary material. Ionic geometries 
were obtained for four different configurations: unligated 
deoxyheme, FeP-d; the heme-CO complex carboxyheme, 
FeP(CO); the heme-02 complex oxyheme, FeP(02); and 
a theoretical planar version of deoxyheme, FeP-p. 

We first discuss the dependence of the iron 3d sub- 
space occupancy rid on the Hund's coupling parameter J 
(Fig. La). We emphasize that the expectation value of 
the occupancy rid of the iron 3d sub-shell is not con- 
strained to integer values in DFT and DFT+DMFT, 
since the iron occupation is a local observable, and hence 
does not commute with the Hamiltonian and is not con- 
served and there are valence fluctuations. 

In the typical region of physically meaningful values of 
the Hund's coupling for iron 3d electrons, J « 0.8 eV, [22] 
we find a very sharp dependence of the electronic density 
on J. In fact, J « 0.8 eV places heme directly in the tran- 



sition region between low-spin states and the = 5 e 
fully-polarized state obtained for large Hund's coupling. 
We note that our results are weakly dependent on the 
choice of the Coulomb repulsion U (see sup. material). 

In Fig. l.b, we show the effective quantum spin num- 
ber, which is associated to the norm of the angular spin 
vector S by the usual relation |S| = y^s(s + 1). The 
spin s shows characteristic plateaux as a function of 
the Hund's coupling at the semi-classically allowed val- 
ues of the magnetization (corresponding to pure doublet, 
triplet, quartet, and quintet states). A fully-polarized 
state is recovered for sufficiently large Hund's coupling, 
as expected. 

At J = 0.8 eV, and almost irrespective of ligation and 
doming, we find that heme has a spin expectation value 
of s « 1.5 corresponding to a quartet state in a semi- 
classical picture. Our results indicate that the true many- 
body wave-function of FeP-d is thus an entangled super- 
position of triplet and quintet states. The proposition 
that heme might be in an entangled state was pointed out 
early [23] in the context of a Pariser-Parr-Pople model 
Hamiltonian, and is confirmed by our DMFT calcula- 
tions. In particular, this accounts for the striking differ- 
ences obtained experimentally for very similar porphyrin 
systems, e.g. it was found that unligated FeP is a triplet 
[24] in the tetraphenylporphine configuration, a triplet 
with different orbital symmetry in the octaethylporphine 
configuration [25], and a quintet in the octamethyltetra- 
benzporphine configuration [26]. The strong dependence 
of the spin state with respect to small modifications in 
the structure is consistent with an entangled spin state. 

In our calculations, we find that both oxyheme and 
carboxyheme adopt a low spin state for J < 0.25eV and 
larger multiplicities in the physical region of J « 0.8eV, 
while in both cases the spin state is very close in character 
to that of unligated deoxyheme. Significantly, we observe 
only subtle differences between FeP(02), FeP(CO) and 
FeP-d for J = 0.8 eV, while the DFT and DFT+£7 treat- 
ment yields ground-states for carboxyheme and oxyheme 
of pure closed-shell and open-shell singlet configurations, 
respectively [6, 7, 9]. 

Moreover, we find that the symmetry of the highest 
occupied molecular orbital (HOMO) of FeP-d, as esti- 
mated from the real-space spectral density of the promi- 
nent feature below the Fermi level, is highly dependent 
on the Hund's coupling J. In particular, for J = eV, 
the HOMO is an admixture of orbital characters (see ver- 
tical labels in Fig. l.a). However, the Hund's coupling 
drives a rather complex orbital selection, such that for 
the region of greatest interest, J ^ 0.8 eV, the HOMO 
predominantly exhibits d 3 ^2_ r 2 symmetry. The orbital 
selection process also induces a pinning of the Fermi den- 
sity to the quantum impurity, such that it is delocalized 
for J = eV (see Fig. l.c), while for J = 0.8 eV (Fig. l.d) 
it is instead localized to the iron 3d sub-shell. 

In our view, this relates to the Fe-O-0 angle obtained 
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FIG. 2: Valence fluctuations: a) Von Neumann entropy 
A obtained by DFT+DMFT for FeP-p (circles) and FeP-d 
(squares). Histograms of the dominant electronic configura- 
tions for FeP-d for b) J = eV and (c) J = 0.8 eV. The pie 
wedge labelled other contains configurations with a weight 
smaller than 3%. The iron 3d spin S z and iron 3d occupancy 
rtd of the dominant configurations is indicated. 

in FeP(0 2 ) [27]. Indeed, the bent geometry of FeP(0 2 ) 
can be explained by a favorable interaction between the 
p*-orbital of the 2 and the d^ z 2_ r 2 -orbital on Fe [27]: 
the O2 p*-orbital is closer in energy to d 3z 2_ r 2 compared 
to the p*-orbitals in CO, and hence it gains more energy 
by bending, which increases the overlap. For FeP(CO) 
the situation is opposite, and there is no stabilization 
gained by bending [27]. On the contrary, the bending 
in FeP(CO) is induced by the strain of the protein and 
it reduces the binding energy. Naively, the orbital se- 
lection of the d 3z 2_ r 2 orbital is hence expected to go 
in the direction of curing the strong O2 and CO imbal- 
ance. Moreover, the charge localization at the Fermi level 
suggests that other artificial binding between the non- 
metallic atomic orbitals of heme and strong electronega- 
tive O2 will not be obtained, and hence will protect heme 
from undesired charge transfer. 

We now discuss the degree of quantum entanglement 
exhibited by FeP-d and FeP-p (see Fig. 2. a). We com- 
puted the von Neumann entropy A = — tr (pd log(p^)), 
where pd is the reduced finite-temperature density- 
matrix of the iron 3d impurity subspace, traced over the 
states of the AIM bath environment. The entropy quan- 
tifies to what extent the wave-function consists of an en- 
tangled superposition. 

We observe that the entropy rises sharply at J ~ 
0.25 eV, corresponding to the transition from the dou- 
blet spin state to the triplet /quintet entangled state. As 
expected, the entropy is small in the low-spin region 
(J < 0.25 eV) and also in the fully-polarized limit. At 
J = eV (Fig. 2.b), we find that the dominant config- 
uration consists of the doublets (d 3z 2_ r 2) 2 (d xy ) 2 (d xz ) 2 , 
with a single electron in the d x 2_ y 2 orbital. The lat- 
ter hybridizes strongly with the nitrogen 2p orbitals, but 
all other orbitals are mostly filled or empty, so this con- 
figuration is, essentially, a classical state with a finite 
magnetic moment. 

At larger J values, however, such as J — 0.8eV 
(Fig. 2.c), all orbitals are partially filled, and an increas- 
ing number of electronic configurations, with different 
valence and spin, contribute to the statistics, and thus 
the iron impurity wave-function is fluctuating. Although 
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FIG. 3: Energetics: a) Difference in CO and O2 binding 
energies AAE. The binding to CO is always favored, however 
the imbalance is strongly reduced for J > 0.5eV. b) Total 
energy of FeP as a function of J. The minimum of the total 
energy is obtained for J = 0.9eV. 

the valence fluctuation are captured to some extent at 
the DFT level (Ajjft ~ 0.75), we find that many body 
effects contribute significantly to the entropy. 

Our results indicate that as FeP-d and FeP-p molecules 
approach a regime with large entanglement for J « 0.5, 
with a concomitant orbital selection close to the Fermi 
level. The orbital selection close to the Fermi level in 
turn induces a charge- localization effect. The latter ef- 
fect of the Hund's coupling can be understood with a 
simple picture: a large Hund's coupling partially empties 
the d 3z 2_ r 2 orbital and brings the weight of this orbital 
closer to the Fermi level, thereby reducing the hybridiza- 
tion between the iron 3d states and the nitrogen 2p states 
close to the Fermi level. The subtle interplay between the 
charge-localization induced by the Hund's coupling (or- 
bital selection close to the Fermi energy) and the dereal- 
ization induced by strong correlations (the tendency for 
electrons to escape the iron 3d orbitals in order to reduce 
the Coulomb energy) is captured by the DFT+DMFT 
methodology but is absent in Kohn-Sham DFT. We em- 
phasize that these ingredients are paramount to an es- 
timation of the charge transfer and binding properties 
between the iron atom and the ligand in oxyheme and 
carboxyheme. 

Let us next discuss the effect of the Hund's coupling 
with respect to the unrealistic imbalance between the 
binding energies of CO and O2 obtained by DFT. The 
binding energy is defined as: AE = E(FeP(X)) — 
(E(FeP) + E(X)), where X=CO or X=0 2 . The differ- 
ence between the binding energies AE(CO) — AE{02) 
is obtained by: AAE = AE C o - A£o 2 • For J = OeV, 
we find that the binding to CO is dramatically favoured, 
when compared to the binding to O2 (Fig. 3. a): the dif- 
ference in binding energies is of the order of 5eV. Al- 
though the binding to CO is favoured for all values of J, 
we find that it is dramatically improved for J > 0.5eV, 
and is reduced down to leV. This suggests that other 
effects might be important to reduce further the CO/O2 
imbalance, such as that the effect of the protein via the 
bending of the Fe-C-0 angle [9]. 

It is also worth noting that we find that the total 
energy of the molecule is minimized for J = 0.9eV 
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FIG. 4: Optical measurements: Optical conductivity of 
FeP-d (bold line) and FeP(CO) (solid) and FeP(0 2 ) (dashed 
line). The vertical arrow indicates the energy of the experi- 
mental peak associated to the Fe-N charge transfer. Inset: ex- 
perimental measurements [28] for unligated heme (bold line), 
oxy- (dashed line) and carboxy- (solid) heme are shown for 
comparison. J = 0.8eV was used for all calculations above. 

(Fig. 3.b), suggesting further that the heme molecule is 
particularly well suited to host metallic d atoms, which 
tend to have a large screened J interaction when hybri- 
dising to light elements such as nitrogen or oxygen. 

We now move to our calculations of the optical absorp- 
tion spectra of heme (Fig. 4). Our theoretical absorption 
spectra, shown in Fig. 4, are in reasonable agreement 
with experimental data [28], in particular for the optical 
transitions at uo ~ 2 eV. We attribute this spectral fea- 
ture to charge-transfer excitations from iron to nitrogen- 
centered orbit als. The spectrum is dominated by the 
characteristic porphyrin Q-bands (those at « 2eV), and 
Soret bands [29] (at « 4eV). Our results offer insight 
into the infrared absorption band present at ~ leV, in 
our calculation, and observed in experiments at 0.6 eV 
[30]. This infrared peak is described, in our calculations, 
as arising from transitions between the d 3z 2_ r 2 spectral 
feature (HOMO) below the Fermi level and the LUMO 
(quasi-degenerate d xz and d^) above the Fermi level. 

Interestingly, we find that the infrared optical weight 
in unligated heme, associated with d-d transitions and 
present in FeP-d, is absent in the planar theoretical 
model FeP-p. Hence, the symmetry breaking associated 
with the doming effect of the iron-intercalated porphyrin 
macrocycle permits d-d optical transitions, and is respon- 
sible for the spectral weight in the infrared regime. We 
note that experimental spectra for FeP(CO) and FeP(02) 
exhibit a double peak structure at uj w 2eV, absent from 
our calculations done at J = OeV, but recovered for 
J > 0.8eV. The best agreement with the experimental 
data is obtained for J = 0.9eV. Finally, we extended our 
calculations to the time dependence of the magnetization 
of the iron atom after an initial quench in polarization 
(see sup. material). We propose that time-resolved spec- 
troscopy may be used as a sensitive probe for the ligation 
state of heme. 

In conclusion, we have carried out linear-scaling first- 
principles calculations, in combination with DMFT, on 



both unligated and ligated heme. We have presented a 
newly-developed methodology applied to a molecule of 
important biological function, exemplifying how subtle 
quantum effects can be captured by our methodology. 
In particular, we have found that the Hund's coupling 
J drives an orbital selection process in unligated heme, 
which enhances the bonding in the out-of-plane direction. 
The von Neumann entropy quantifying valence fluctua- 
tions in the iron 3d subspace is large for the physical 
values of J « 0.8eV. This scenario sheds some light on 
the strong CO and O2 binding imbalance problem ob- 
tained by extracting the binding energies in simpler zero 
temperature and J = OeV DFT calculations. The dif- 
ference in binding energies is dramatically reduced for 
physical value of J ~ 0.8. The smaller remaining im- 
balance might be further explained by the strain energy 
contained in the protein structure [9] or by the contri- 
bution from the entropic term. Finally, the relevance of 
a finite Hund's coupling in heme is confirmed by the to- 
tal energy extracted from the DFT+DMFT of unligated 
heme, which shows a minima for J = 0.9eV. 

We have proposed a new mechanism for ligand binding 
to heme based on an orbital selective process, on this ba- 
sis, a scenario which we term bonding determined by local 
valence fluctuations. Finally, we have obtained a reason- 
able agreement between experimental and our theoreti- 
cal optical absorption spectra, our description accounting 
for the observation of optical transitions in the infrared 
regime and the double peaked structure of the optical 
response at uj ~ 2eV. 

At the time of writing, we became aware of related 
application of DMFT to an organometallic crystal [31]. 
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In this supplemental material, we study the dynamical magnetic response of heme, in various 
ligation and doming configurations, by means of out-of-equilibrium calculations performed by per- 
turbing the DFT+DMFT ground state wave-function and propagating in the Keldysh time domain. 
Based on the rather strong ligand specificity observed in these simulations, we propose a potentially 
sensitive method for profiling of heme binding, using two-photon circularly-polarized spectroscopy 
Here, we also describe the theoretical methods used in the work presented in our Letter in detail, and 
justify the range of parameters used in our study, i.e. we explore the dependence of our calculations 
on the choice of the Coulomb repulsion and the temperature. A simple toy model for calculating 
the entropy is proposed, and compared to the DMFT calculations, in order to illustrate the impor- 
tance of many body effects. We furthermore justify the importance of the Hund's coupling based 
on phenomenological grounds, by carryout a comparison between the theoretical optics at various 
values of the Hund's coupling with experiments. Out-of-equilibrium calculations are also discussed. 



In this work, we have carried out a detailed theo- 
retical study of the electronic structure of the heme 
molecule by means of a combination [1, 2] of linear-scaling 
density-functional theory (DFT) with the dynamical 
mean-field theory approximation (DFT+DMFT) [3, 4], a 
model which includes local dynamical, finite-temperature 
and multi-determinental effects, for given Hamiltonian 
parameters, at an effectively exact level. Although 
DFT+DMFT has been primarily applied to solids to 
date, here we apply our real-space DFT+DMFT imple- 
mentation to a substantial molecular system, extending 
the scope of DFT+DMFT applicability to biology for the 
first time to our knowledge. We performed a study of the 
impact of varying the Hund's exchange coupling strength 
J. Although the value of J in iron based molecules and 
solids is expected to be J ~ 0.8 eV, [5, 6] we are interested 
in this work to study what are the physical ingredients 
introduced by J as it is increased from zero. In particu- 
lar, it was suggested that the Hund's coupling might fix 
the strength of the correlations [7] , or act to increase the 
total magnetic moments within the localized effective im- 
purity subspaces treated using DMFT and, in doing so, 
break degeneracies between the localized orbitals. 



The effects on the electronic structure due to the di- 
atomic axial ligands O2 and CO are considered, and test- 
ing is further carried out with respect to changes in the 
Hubbard U parameter and temperature. Our study fur- 
thermore includes two experimental signatures of ligand 
binding, namely the optical conductivity discussed in our 
Letter, and the transient magnetic response with which 
we discuss at the end of this supplemental material. 



METHODOLOGY: DENSITY-FUNCTIONAL 
THEORY CALCULATIONS USING OPTIMIZED 
WANNIER FUNCTIONS 

Kohn-Sham density-functional theory (DFT) [8, 9] was 
used in this study to provide a reasonable, conveniently- 
computed zero-temperature, single-determinental start- 
ing point for the DFT+DMFT self-consistency proce- 
dure, within the localized subspace treated for many- 
body effects using DMFT, as well as a sufficiently de- 
tailed description of the delocalized electron "bath" com- 
prising the remainder of the system. Relatively simple 
local and semi-local approximate functionals for DFT 
often provide an acceptably refined description of the 
ground- state electronic structure, and even the quasi- 
particle band-structure, associated with orbitals of large 
spread or no magnetic order. This observation forms the 
basis of the DFT+DMFT approximation, in which the 
dominant computational effort is instead focused on solv- 
ing for the many-body wave-function constructed from 
the the remaining, localized, possibly spin-polarized or- 
bitals adjacent to the Fermi-level. 

The ONETEP linear-scaling DFT and DFT+Z7 
code [10-12] was used to solve for the DFT ground- state 
subspace Green's function in this work. The ONETEP 
method is based on direct minimization of the total- 
energy with respect to the single-particle density-matrix, 
and is particularly advanced in terms of its accuracy 
equivalent to that of a plane-wave method, which is ar- 
rived at by means of an in situ variational optimiza- 
tion of the expansion coefficients of a minimal set of 
spatially-truncated Nonorthogonal Generalized Wannier 
Functions [13] (NGWFs) (with non-trivial overlap 

Sap — (<t>a 1 0/?))- The basis for this expansion is a set 
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of systematically improvable Fourier- Lagrange, or psinc 
functions [14], which is equivalent to a truncated set 
of plane- waves. We express the single-particle density- 
matrix in the separable form [15, 16] given by 

po(r,r') = Y,<Pa(r)K^4> (v'), (1) 

where K a ^ = (cj) a \po\(j)^) is known as the density- kernel, 
the tensor representation of the single-particle density 
operator po [58]. The contravariant NGWF duals {<j) a } 
are defined as those which satisfy the biorthonormality 
relation, where the overlap matrix acts as metric tensor, 
defined by 

(<f> a \<f>p)=6%, with \4> a ) = {S- 1 ) a0 \^)- ( 2 ) 

The use of a minimal, optimized Wannier function rep- 
resentation of the density-matrix allows for the DFT 
ground state to be solved with relative ease in large sys- 
tems, particularly in molecules where their explicit trun- 
cation implies that the addition of vacuum does not in- 
crease the computational cost. Furthermore, it allows 
for the full Kohn-Sham Green's function to be computed 
via Hamiltonian inversion with a manageable computa- 
tional overhead. Thus, we do not make any projection 
over bands or orbitals prior to generating the full Green's 
function, and only project after inversion when comput- 
ing the matrix elements of the subspace Green's function. 



place of the core electrons for all ions. The semi-core 
states of iron were included in the valence, so that the 
(3s,3p,3d,4s,4p) sub-shells were were retained for explicit 
treatment, obviating the use of a non-linear core correc- 
tion. For the lighter element, (2s, 2p) were retained in the 
valence for C, N, and O, along with (Is) for H. Hence, we 
used 13 variationally-optimized NGWFs to describe the 
Fe electrons, 4 NGWFs for each C, N, and O atom, and 
1 for each H atom, giving a total of 213 NGWFs (spin- 
degenerate orbital pairs) for the unligated heme system 
and and 221 each for the two ligated systems. 

Ionic geometries were obtained from the Protein Data 
Bank: unligated deoxyheme (dome-shape, FeP-d) was 
extracted from a human deoxyhaemoglobin structure ob- 
tained by x-ray crystallography (PDB key 2HHB [20]), 
and the missing hydrogens were reintroduced using the 
Jmol package [21]. The Carboxyheme (heme-CO com- 
plex, FeP(CO)) and the oxyheme (heme-02 complex, 
FeP(02)) structures were extracted using the same pro- 
cedure (PDB key 1GCW [22] and 1HHO [23], respec- 
tively). A theoretical, planar deoxyheme (FeP-p) (PDB 
key HEM) is also considered in our calculations for com- 
parison with the experimentally-observed dome-shaped 
structure in order to test how the conformal change be- 
tween the planar and dome-shape affects the electronic 
structure. In the present work, we limit our calculations 
to a single heme molecule, omitting residues which neigh- 
bour it in the haemoglobin protein. 



METHODOLOGY: IONIC GEOMETRIES, 
PSEUDOPOTENTIALS AND VARIATIONAL 
PARAMETERS 

We performed ground-state DFT calculations on heme 
in vacuo by iteratively minimizing the energy functional 
with respect to both the density kernel and the NGWFs, 
using the algorithm described in Ref. [17]. The total- 
energy was converged to within 1 meV per atom, with 
respect to the plane- wave energy cutoff (at a minimum 
of 1300 eV for all reciprocal lattice vectors in our cal- 
culations) and the NGWF cutoff radii (at 6.6 A for all 
species) , and no additional restrictions on the variational 
freedom, such as the density kernel truncation, were in- 
voked. The simulation cell volume was increased until 
(at 4.8 x 10 4 A 3 ) periodic-image interaction errors also 
fell within this total-energy tolerance. 

The DFT calculations were performed, and pseudopo- 
tentials generated, using the PBE [18] generalized gra- 
dient exchange-correlation functional. In line with com- 
mon practice in DFT+DMFT calculations, in which all 
magnetic order is assumed to be confined to the local- 
ized impurity subspaces, spin-degenerate DFT ground- 
state densities were used. Scalar scalar-relativistically 
corrected norm-conserving PBE [18] pseudopotentials, 
generated with the OPIUM package [19], were used in 



METHODOLOGY: DEFINITION OF THE 
LOW-ENERGY, CORRELATED MODEL 

We extended our DFT calculations using the 
DFT+DMFT method [3, 4] in order to refine the descrip- 
tion of strong correlation effects arising due to partial de- 
generacies between the localized orbitals making up the 
iron 3d sub-shell at the core of heme. DMFT allows for 
quantum and thermal fluctuations, absent from the zero- 
temperature, single-determinant Kohn-Sham DFT calcu- 
lations, to be accurately described. The heme molecule 
was represented within DMFT by an Anderson impurity 
Hamiltonian (AIM) [24] , in its Slater-Kanamori form, the 
ground-state of which was solved for by using a finite- 
temperature Lanczos solver [25]. Since only a single im- 
purity site (3d orbital subspace) is present, the system 
becomes crystal momentum independent in the molecu- 
lar limit, and since the Kohn-Sham Green's function is 
inverted in full before projection onto the impurity sub- 
space, the Anderson impurity mapping is effectively ex- 
act, so that the DFT+DMFT algorithm is almost entirely 
converged in a single Green's function self-consistency 
step. We note that a Kondo model was used early in the 
literature to study heme, along similar lines than ours, 
by means of a model Hamiltonian [24]. In our study, the 
AIM model is determined and obtained by first-principle 
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calculations. 

A spherically symmetric trial impurity subspace was 
defined by a spanning set of iron 3d orbitals, the or- 
thonormal set {ip m } produced by solving the spherically- 
symmetric Schrodinger equation subject to the norm- 
conserving iron pseudopotential with an appropriate con- 
fining potential, the within a truncation sphere of radius 
6.6 A. The same procedure was used to generate the 
initial guesses for the NGWFs during the initialization 
of the DFT calculation, so that the trial impurity sub- 
space formed a proper subspace of the initial guess for 
the Kohn-Sham Hilbert space. Since the latter is op- 
timized as the energy is minimized with respect to the 
NGWFs, so, at convergence of the DFT algorithm, the 
Hubbard projectors finally chosen to span the impurity 
subspace of DFT+DMFT were the so-called symmetry- 
adapted Nonorthogonal Generalized Wannier Functions 
(SNGWFs), defined by 

\<Pm) = |0a)(0 a bm) = |0a> (M<Pm)- ( 3 ) 

Thus, the final impurity subspace is a proper subspace 
of the converged Kohn-Sham Hilbert space, a necessary 
condition for a strictly causal self-energy, but retains the 
3d symmetry of the numerical atomic orbitals (the gen- 
eralization to nonorthogonal projector functions in this 
context is discussed extensively in Ref [26]). 

Once the fully- converged DFT energy minimization of 
was carried out, for each heme strain, the full Green's 
function was initially computed in the finite-temperature 
Matsubara representation. Noting that the inverse of a 
doubly- covariant tensor is a doubly-contravariant tensor, 
the full Green's function is generally expressed and com- 
puted, in terms of the the Kohn-Sham Hamiltonian H, 
as 

G al3 (iw„) = ((iw„ + v)S a p - H aP - E^)- 1 . (4) 

Here, \i is the chemical potential, set to the average of 
the highest occupied and lowest unoccupied Kohn-Sham 
orbital energies, and X is the self-energy tensor gener- 
ated by the DMFT algorithm, where XI = in the first 
instance so that the initial full Green's function is the 
Kohn-Sham Green's function Go- Green's function sam- 
pling at 400 Matsubara frequencies provided adequate 
convergence of the properties of interest. We performed 
this matrix inversion, as well as all matrix multiplications 
involved in the DMFT algorithm, on graphical compu- 
tational units (GPUs) using a tailor-made parallel im- 
plementation of the LU decomposition using the CUDA 
programming language. This provided a crucial improve- 
ment in the computational feasabililty of our calculations. 
Following inversion to find the the full Green's function, 
the Kohn-Sham subspace Green's function Go is given 
by its projection onto the impurity subspace, where it 



has the matrix representation 

Gomm' fan) = (<Pm\G (i(jJ n ) \<fm') (5) 

where m and m' run over the five iron 3d SNGWF projec- 
tor functions (in real cubic-harmonic notation: d x 2_ y 2 1 
ds z 2_ r 2, d yzi d XZl d xy ), a and j3 are the indices for 
the NGWFs, and the matrices NGWF-projector over- 
lap matrices are defined as Vain = {(t>a\^m) and Wmi = 

(<Pm\<l>a)- 

In practice, in order to imbue the SNGWF Hubbard 
projectors with a more plausible physical interpretation, 
a real-space rotation of the functions was carried out [27] 
in order to better align their lobes. The subspace pro- 
jected Green's function is thus transformed to 

G rot = U+GU, (6) 

where U is the 5x5 rotation matrix in cubic harmonic 
space, corresponding to a rotation in R( 3 \ and G rot 
is the rotated subspace Green's function passed to the 
DMFT solver. The rotation matrix U is chosen, prag- 
matically, such that the e x and e y axes are those, in an 
averaged sense, which point towards the four nitrogen 
atoms surrounding the iron-centered impurity subspace, 
with e z directed out of the porphyrin plane. 

Electronic correlation effects beyond the capacity of 
the approximate exchange-correlation functional, those 
arising due to interactions between particles within the 
impurity subspace and finite- temperature effects, are 
explicitly described in our DMFT calculations by the 
Slater-Kanamori form of the Anderson impurity Hamil- 
tonian [28, 29], specifically 

7-Lu = U y^n mt n m ^ + (U f - ^) ^ n m n m / (7) 

m m>m' 

-j y, (2s n ,s m / + (4, t 4 4 d m / t d m ^)). 

m>m / 

In this, the first term describes the effect of intra-orbital 
Coulomb repulsion, parametrised by U, and the second 
term describes the inter-orbital repulsion, proportional to 
U' , which is renormalized by the Hund's exchange cou- 
pling parameter J in order to ensure a fully rotationally 
invariant Hamiltonian (for further information on this 
topic, we refer the reader to Ref. [30]). The third term is 
the Hund's rule exchange coupling, described by a spin 
exchange coupling of amplitude J. S m denotes the spin 
corresponding to orbital m, so that S m = |dj ns tr ss /d ms /, 
where a is the vector of Pauli matrices indexed by s and 
s' . We note that the Slater-Kanamori form of the ver- 
tex interaction is especially well suited to capture the 
multiplet properties of the low energy states [30]. The 
Slater-Koster [31] interaction is another alternative, but 
is mostly used to describe solids, and is not well suited 
to the case of a molecule. 
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In this this work, unless where otherwise stated, we 
used U = 4 eV for the screened Coulomb interaction [32] 
and we explored the dependence of several observables 
on the Hunds coupling (in the range J = — 2.5 eV). 
Our DMFT calculations were carried out at room tem- 
perature, T = 294 K, again except where otherwise 
stated. The dependence of our results with respect to 
the Coulomb repulsion U and the temperature T is also 
discussed in the next sections. 



METHODOLOGY: DYNAMICAL MEAN-FIELD 
THEORY TREATMENT 

The Hamiltonian (7), in combination with the expres- 
sions (4) and (5), defines the many-body impurity prob- 
lem that we solve using the DMFT algorithm [3], which 
updates the impurity Green's function G. The DMFT 
establishes a self consistent mapping from the correlated 
atoms of the initial solid or molecule and a smaller local 
problem, the Anderson Impurity Model (AIM), which 
is used to obtain the self energy within this projected 
subspace [3]. The mapping is carried out self consis- 
tently, and the obtained local Green's function of the 
AIM converges to the Green's function of the larger space 
projected onto the correlated atom. At the level of the 
AIM, this model describes an impurity connected to a 
bath by the so-called hybridization function. The bath 
of the AIM models the surrounding environment of the 
correlated atom in the solid or the molecule, and the hy- 
bridization function describes how electrons are dynam- 
ically transferred from and to the bath to the impurity. 
Hence, at the AIM level, the uncorrelated part of the 
Heme molecule is described by the bath, and the hy- 
bridization between the d orbitals of the iron atom and 
the N atoms is described by the hybridization function 
in the AIM. 

We define only a single impurity subspace in calcula- 
tions on heme, since there is only one transition-metal 
ion present, and so the impurity self-energy 

E = Go 1 - G 1 (8) 

resulting from the DMFT algorithm is said to be exactly 
local, in the sense that there are no pairs of impurity sites 
for which to consider site off-diagonal self-energy matrix 
elements. In this particular case, there is no feedback 
from the self energy to the hybridization function, and 
the DMFT converges after one iteration. It is hence not 
a mean-field approximation in this particular case, but 
turns out to be exact. However, the methodology de- 
scribed in this work can be applied to molecules with 
many correlated ions with no further modifications. 

The hybridization matrix A within the AIM impurity 
subspace is given, formally, by 

A(iou n ) = (iu n + fi) 6 - 53 - E imp - GT\ (9) 



where here [33-35], the metric tensor on the SNGWFs is 
given by 



6 = (ws _1 v) 1 . 



(10) 



This metric is in general non-trivial, i.e., O ^ 1 and 
so the SNGWFs are nonorthogonal and not identical to 
their duals, even if their parent atomic orbitals form an 
orthonormal set, if the trial impurity subspace does not 
form a proper subspace of the converged Kohn-Sham 
Hilbert space. However, for this particular case of study, 
we found however that the deviation of O from the iden- 
tity matrix was small (within 0.1%). 

The high-frequency part of the hybridization function, 
jnmp _ A(io; n — » oo), defines the impurity energy levels 
in the zero- hybridization "atomic" limit, so that the hy- 
bridization matrix A(icj n ) exhibits the correct physical 
decay proportional to l/iu n . We tested that the implied 



G- 1 ^)" /u> 



holds, 



limiting condition O = lim^- 
up to a high precision, ensuring that the self-energy is 
physically meaningful. It can also be straightforwardly 
obtained by doing a high frequency expansion of the 
Green's function that: 



E im P = OW (S^HS -1 ) vo. 



(11) 



The self-energy T, is thus obtained by solving the 
Anderson impurity model (AIM) defined by the hy- 
bridization (9) and the interaction Hamiltonian (7), using 
a finite-temperature Lanczos DMFT algorithm [36-38]. 
The Lanczos solver uses a finite discretization of the hy- 
bridization (9), suffering finite size effects, yet allows for 
the orbital off-diagonal elements of the hybridization (9) 
to be considered on an equal footing to the diagonal ele- 
ments. Crucially, this mitigates some of the dependence 
on the spatial form and orientation of the projector func- 
tions (SNGWFs) used to define the impurity subspace. 

Having solved for impurity self-energy given by the 
AIM, X) rot , it is rotated back to the original system of 
coordinates, to give 



(12) 



We then "up-fold" it to the Kohn-Sham Hilbert space, 
using its separable form in terms of NGWFs, that is 



^a(3 — V a 



v dc o r 



w, 



m'p- 



(13) 



A separable form of the up-folded self-energy enforces its 
causality, as shown recently in Ref [39], (the local self- 
energy for the impurity subspace being causal by con- 
struction, so that E mm /(zcj n ) < 0, for all ioj n ,m,m . 
We carefully verified that our computed self-energies were 
causal, a necessary condition for positive definite spectral 
functions and the physicality of computed observables 



5 



generally. In this work, we used the canonical form of 
the double-counting potential Vd c , given by 



Vdc = u a 



J 



(n d - 1) , 



(14) 



where rid is the invariant total occupancy of the impurity 
subspace, 



2 

rid = ~ I duj G m7n > (uj) O r ' 
Jo 



7T 



(15) 



Here, the parameter U av is the averaged repulsion, re- 
lated to the intra-orbital and inter-orbital repulsion, and 
computed using the expression [40] 



jja 



U + 2 (7V deg - 1) u> 
2N deg - 1 



(16) 



where N deg = 5 is the number of localised orbitals (SNG- 
WFs) spanning the impurity subspace. The double- 
counting correction is intended to subtract the interac- 
tions within in the impurity subspace that are already 
included, to some extent, in the Kohn-Sham Hamilto- 
nian H. This is carried out in an approximate, averaged 
fashion, since the exact form of double-counting correc- 
tion required for a given density, exchange-correlation 
functional, and set of projecting orbitals, is not known. 
A good test of the validity of the double-counting cor- 
rection, verified in our calculations, is the independence 
of the impurity chemical potential with respect to the 
Coulomb interaction parameter U. The equations (4), 
(5), (9), and (13) together form the Green's function 
self-consistency cycle used in our DFT+DMFT calcu- 
lations, which is iterated until convergence is obtained. 
In the calculations presented in this work, for reasons 
aforementioned, the mapping onto the AIM is effectively 
exact, and the DFT+DMFT self-consistency cycle con- 
verges rapidly, being already close stationary after a sin- 
gle iteration. 



METHODOLOGY: IRON DENSITY AND 
MAGNETIZATION, ENTANGLEMENT AND 
ENTROPY 

Once the AIM problem is solved with the finite temper- 
ature Lanczos algorithm, we have access to the ground 
state and the low energy excited states. The eigenstates 
of the many-body Hamiltonian are written in the Fock 
space of the impurity connected to the bath. The basis 
states of the Fock space span all the possible electronic 
configurations, hence an eigenstate of the hamiltonian is 
a superposition of all possible quantum states of the im- 
purity and the bath. It is important to note that other 
simpler approaches, such as the so-called Configuration 
interaction [41] only consider a very restricted basis set 
where a few electrons are excited from the impurity to 



the bath. In our approach, this is not the case, and all 
possible configurations are taken into account. 

Observables related to the impurity model, such as the 
impurity density, can be obtained by the finite temper- 
ature thermal average. The trace of an operator A is 
defined as: 



52e-P Ei (i\A\i) 



(A) 



(17) 



where Z is the partition function, /3 = 1/T is the inverse 
temperature, and the index i runs over the low energy 
states. The impurity density operator nd is simply de- 
fined as: 



n d 



E4 



(18) 



where c\ a (c ia ) creates (destroys) an electron in the or- 
bital i of the iron d manifold with spin a. The thermal 
average of n d is given by: 



E e-^(i\ci a c acr \i) 

i \ i,ac,cr 

(nd) = ^ 



(19) 



where a denotes the 5 d orbitals, i runs over the eigen- 
states of the AIM, which is determined self-consistent ly 
by the DMFT, E{ is the eigen energy of the quantum 
eigen state \i). We note that the summation over the 
eigenstates \i) involves states with different number of 
electrons, and hence is not an integer number, but can 
take any rational values. Since upon DMFT convergence 
the impurity Green's function converges to the Green's 
function of original problem, the density nd of the im- 
purity model corresponds to the electronic occupation of 
the iron atom in heme. 

We note that, although the total number of electrons, 
and global observables more generally, are associated 
with operators that commute with the Hamiltonian and, 
hence, are conserved quantities, the density at the iron 
site is a local observable and can fluctuate strongly. 

In particular, we find for J = eV that DFT+DMFT 
yields an impurity occupancy of 6.8 e for deoxyheme, 
(both for FeP-p and FeP-d), with larger subspace oc- 
cupancies for ligated hemes (n^ = 7.0 e for FeP(CO) 
and nd = 6.9 e for FeP(02)). This is in good agreement 
with earlier DFT+U = 6 eV calculations [42], which gave 
rid = 6.6 e for deoxyheme and rid — 6.8 e for carboxy- 
heme. We note that doming of the heme molecule has 
no effect on the subspace occupancy at small Hund's cou- 
pling, although a significant dependence of the occupancy 
on the doming is activated beyond J ~ 0.3 eV. 

In our calculations we found that the DMFT solution 
remains paramagnetic, although it was allowed for the 
possibility to spontaneously form a magnetic moment 
(spin symmetry broken state). However, the low energy 
states are in a quantum superposition of polarized states, 
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giving a fluctuating magnetic moment at the iron site, 
which can be quantified as [43]: 

S = yJ(SS)-((S))* = (20) 

Where S is the obtained magnetic moment of the iron 
atom, S the total spin operator at the iron site. The spin 
s is obtained as usual by the relation: 

5 = s(s + l). (21) 

We note that the eigenstates obtained by the Lanczos 
are written in the Fock space which contains both the 
impurity orbitals and the bath orbitals, and although the 
iron density fid operator is defined in terms of d orbitals, 
the low energy states are defined in the basis of both 
the impurity and the bath orbitals. The entanglement 
between the impurity and the bath can be quantified by 
the reduced density matrix of the impurity p, defined as: 

p = J2e-^Tr B \i)(i\, (22) 

i 

where Tr# denotes the partial trace over the bath de- 
grees of freedom, and \i) are the low energy states of the 
problem. 

In particular, the eigenvectors v& of p provide a sim- 
plified description of the many body states spanned by 
the impurity: they give a cartoon picture in terms of 
electronic orbitals for the impurity at low energy. The 

eigenvalues of p, are normalized weights (J2^k — 1) 

k 

which give the probability to measure the corresponding 
many body state. We note that: 

• each v/e is a superposition of many body states of 
the Fock space, so it is a quantum superposition of 
electronic occupations of the impurity orbitals 

• there is in general not only one dominant A&, but 
a collection of eigenvalues with similar amplitudes. 

The Von Neuman entropy, defined as: 

A = -^5> fc ln(A fe ) (23) 

k 

measures the entanglement between the impurity and the 
bath. When a number of with similar amplitudes 
starts to proliferate, the entropy becomes large. In the 
particular case when Ai = 1, the impurity is in a pure 
state and the associated entropy A = 0. In the latter 
case, the quantum state of the impurity is well defined, 
and a spin and a charge can be used to label the ground 
state of impurity problem. 

As discussed in the last section, since in the case of 
DMFT applied to heme there is only one correlated site, 
the impurity observables, such as the impurity density 



rid or the entropy A, readily provide respectively the iron 
charge density and the entropy associated with the iron 
atom in the heme molecule. Hence, if the impurity in the 
AIM is in a pure state (A = 0), one might associate to 
the heme molecule a charge density and a spin. For the 
general case A > 0, the iron atom in the heme molecule 
does not have a fixed spin and valence. This is essentially 
related to the fact that the iron density and spin are not 
global observable, and hence they can fluctuate strongly. 

We emphasize that the valence fluctuations are not 
only captured by DMFT, in particular local observables 
also fluctuate in density functional theory. However the 
scheme above, only valid for a many body description 
(Fock space), does not apply at the DFT level and a dif- 
ferent definition of the entropy needs is called for. To 
characterise the entropy in DFT, the so-called linear en- 
tropy definition needs to be used [44] . The linear entropy 
of the reduced density matrix p re d, 

L = Tr(p red )-Tr(p 2 red ) (24) 

may be used as a measure of entanglement for a pure 
state. Since the reduced density matrix is normalized 
(Trp re d = 1), we obtain: 

L = 1 - Tr (p 2 red ) (25) 

With this definition of the entropy, suitable for DFT 
calculations, we found that L = 0.757 in FeP-p and 
L = 0.765 in Fep-d. The entropy at the DFT level 
is hence much smaller than the DFT+DMFT entropy, 
which is in the range A = 2.5 — 4.2. This is explained by 
the absence of the many-body excitations induced by the 
correlations (U and J) which contribute significantly to 
the entropy. 

We note that the linear entropy is an approximation of 
the Von Neumann entropy, which is often used in DFT 
calculations since it does not require a direct inversion of 
the full density matrix p. Although in principle the DFT 
entropy obtained for the exact density matrix p would 
capture the effect of quantum entanglement at zero tem- 
perature, in practice approximation are made for the ex- 
change energy functional E^, which in general do not 
fully capture quantum fluctuations (the extent to which 
DFT reproduces the entanglement can be studied in sim- 
ple problems where the exchange functional is known ex- 
actly, such as in the Hooke's atom [44]). 

A simple many body version can be constructed from 
the DFT reduced density matrix p re d in the spirit of the 
configuration interaction approach, where we have many 
body states. In the reduced subspace, if we neglect strong 
off-diagonal components in the DFT reduced density ma- 
trix p r ed, the statistical operator is diagonal, and we can 
associate to each orbital a many body configuration with 
a statistical weight based on the occupation density of 
this orbital n and on the averaged double occupancy of 
this orbital d. In particular, we associate to each orbital 
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[45]: 



Pa = Po|0)(0| + ^2piW)(a\ +pa|2)(2| (26) 



where |0),|a),|2) denotes respectively no electron, an elec- 
tron with spin a, and two electrons in the orbital a. A 
natural estimation of po,Pi and ^2 is: 

• Po = 1 — n + d 

• Pi = n/2 — d 

• P2 = d 

Because of its diagonal form the local von Neumann en- 
tropy is given by A = — ^p a mp a , with a = 0, *f , ^, 2 and 

a 

the total entropy for the five orbitals a can be obtained 
similarly. To obtain an upper bound on the entropy, we 
can consider the case without correlations U = and as- 
sume d = n 2 /4. We find A = 4.72 for FeP-p and A = 4.76 
for FeP-d. These values are of the same order at the max- 
imum of the entropy obtained with DMFT in the regime 
J « 0.8eV. 

In particular, in the simple argument above, we as- 
sumed that the probabilities associated with every or- 
bital a were independent (no correlations). This simple 
picture breaks in the presence of a finite U. Indeed, we 
find that for J = OeV and U = 4eV the entropy is much 
smaller. A large repulsion quite generally reduces the 
entropy and orders the charge [45]. 

In conclusion, to take into account many body effects 
in the presence of correlations in the iron atom (U and 
J), a consistent method is required. DMFT provides a 
framework to carry out these calculations without any 
add-hoc assumptions. 



METHODOLOGY: SPECTRAL FUNCTIONS 
RESOLVED IN REAL-SPACE AND ENERGY, 
THE OPTICAL CONDUCTIVITY AND THE 
TOTAL ENERGY 

The Lanczos DMFT solver used in our study [46] read- 
ily provides information on the real frequency axis. For 
example, the NGWF-resolved spectral density p a p (u) is 
defined, via the retarded Green's function obtained at 
DMFT self-consistency, G, by 



a/3 



G ali (w) - (u) 



2i7T 



(27) 



We note that the Green's function has poles on the real 
axis, such that it is in principle not possible to define 
G(uj). Hence, the Green's function is defined with a small 
but finite broadening factor in, and G(u + irj) is used in 
the formula above. The small shift n plays the role of a 
broadening factors, similarly to the broadening used in 



DFT to compute the density of states at zero tempera- 
ture. 

The DFT+DMFT density of states is then given by 
the imaginary part of the many-body density matrix, so 
that 



p(u>) = $S[p a P (u>)]Sp a . 



(28) 



The real-space representation of the Fermi density, the 
spectral density at the Fermi level is given by 

n F (r) = a (r) 9 [p^ (u> -+ p)] (r) . (29) 

The theoretical optical absorption was obtained in 
DFT+DMFT within the linear-response regime (Kubo 
formalism), including the effect of non-local pseudopo- 
tentials on the velocity operator matrix elements [47] , in 
the no-vertex-corrections approximation [48], where it is 
given by: 



a(u) 



I 



duo 



2-Ke 2 h 

X (p a ? (w'-Lj)vp yP 1 S (0/)v 4a ), 



(30) 



and the factor of two accounts for spin- degeneracy, ft is 
the simulation-cell volume, e is the electron charge, h 
is the reduced Planck constant, / (uj) is the Fermi-Dirac 
distribution, and the many-body density matrix p a ^ is 
that defined by equation (27). We note that the opti- 
cal conductivity will be smoothed out by the broadening 
parameter irj used to obtain the Green's function on the 
real axis. In practice, in is chosen to be small enough 
such that the important features in the optics can be 
identified. 

The matrix elements of the velocity operator along one 
of the direction e x ,e y ,e z , denoted as (j = (x,y,z)) 
are given by 



T J — 

/ a(3 - 



Vnl,T 



\M- (31) 



This expression is general to the nonorthogonal NGWF 
representation [49], used in this work, where the con- 
tribution to the non-interacting Hamiltonian due to 
the non-local part of the norm-conserving pseudopoten- 
tials, [50, 51] represented by V n i, is included. It is worth 
noting that we do not invoke the Peierls substitution [48] . 

Since optical experiments on heme are typically carried 
out in liquid or gas phases, therefore the isotropic part 
of the optical conductivity tensor: 



o~{u) = ^cFa (cj)/3, 



(32) 



is of primary interest. Finally, we note that optical tran- 
sitions between hybridised d orbitals (d-d transition type) 
are forbidden by optical selection rules if the hybridiza- 
tion of the d and p-orbitals happens in a perfect octahe- 
dral symmetry , such as in FeP-p. 



8 



Once the self energy matrix is obtained, it can be used 
to correct the DFT total energy with the DMFT cor- 
rection. The total energy is estimated by the functional 
[52, 53]: 



E = E LDA (p(r)) - + Tr 



HldaG + (Hu) — E D c 
(33) 

where Hjj indicates the many body interaction vertex 
of the DMFT, and the primed sum is over the occupied 
states. The symbol "Tr" indicates the one-electron trace 
for a generic representation and the sum over the Mat- 
subara frequencies iuj n for finite-temperature many-body 
formalism. 

We notice that the total energy within the 
LDA+DMFT scheme is not simply the expectation value 
of this Hamiltonian but it consists of several terms, in 
analogy to the expressions of the usual DFT. The first 
term Eld a contains four different contributions, namely, 
the ones due to the external potential, the Hartree po- 
tential, the exchange-correlation potential, and the sum 
of the Kohn-Sham eigenvalues. However in the spec- 
tral density- functional theory, the Kohn-Sham eigenval- 
ues should be reoccupied with respect to the descrip- 
tion given by the total Greens function. Then we should 
remove the bare Kohn-Sham eigenvalues sum (second 



H 



LDA 



G 



(third term). 



term) and substitute it with Tr 

The interaction term Hjj can be obtained with the the 
Galitskii-Migdal formula [54]: 



(Hu) 



1 



-Tr 



EG 



(34) 



Where fi (G) is the self-energy (Green's function) matrix 
in the NGWFs representation. We note that both the 
self energy and the Green's function are slowly decaying 
function, hence the trace over matsubara frequencies has 
to be done with care. In this work, we followed the steps 
of Refs. [52, 53]. 

Finally, the double-counting correction E^c must be 
introduced, since the contribution of interactions between 
the correlated orbitals to the total energy is already par- 
tially included in the exchange-correlation potential de- 
rived from the DFT. The most commonly used form of 
the double-counting term is [40]: 



E DC = ^N d (N d - 1) 



J 



J2 N dA^da-l) (35) 



where f7 av is the averaged Coulomb repulsion defined in 
equation (16). 

VALIDATION OF THE METHODOLOGY 

In order to validate our approach, we have studied the 
dependence of our results with the other parameters of 



the theory, e.g. the temperature T and the Coulomb 
repulsion U. 

We show in Fig. l.a the iron occupation as a func- 
tion of J for different values of U = 3.5eV, U = 4eV, 
U = 4.5eV, and U = 5eV for FeP-d. We find that the 
results are qualitatively similar, for a 1.5eV variation of 
U (3.5 — 5eV), although the critical value J associated 
with the sharp drop of the iron density is slightly shifted 
from J = 0.8eV to J = leV. We note however that for 
J = OeV, the density is only very marginally affected as 
a function of U. The latter suggest that strong modifi- 
cations of the physical properties are induced principally 
by the Hund's coupling. In our work, we chose U = 4eV, 
in agreement with earlier DFT studies [32] , and we show 
that although variations in U might give some quantita- 
tive modifications at J=0, there is a dramatic change in 
the occupation density at finite J, neglected so far in the 
DFT studies of heme. We also note that it is physically 
better to take into account the Hund's coupling when 
dealing with metallic atoms such as iron, known to have 
a large Hund's coupling. 

We also studied how our results depend on the temper- 
ature T (see Fig. l.b). We compared the data obtained 
at T = 100K and T = 293K. Despite some minor quan- 
titative differences, the sharp transition in the iron den- 
sity is obtained for the same critical Hund's coupling J, 
showing that our predictions are valid for a broad range 
of physically relevant temperatures. This confirms also 
that quantum fluctuations are mostly dominating in this 
regime of temperature, although we cannot exclude that 
for T < 100K we might observe some larger changes. 

In conclusions, we believe that our predictions are not 
dependent on the details of our model, and we show 
that the parameter which has a dramatic impact on the 
physics in the Hund's coupling. 

To support further the justification that the Hund's 
coupling is important for heme, we now turn to the 
discussion of the optics for FeP-d (Fig. 2. a), FeP(CO) 
(Fig. 2.b) and FeP(0 2 ) (Fig. 2.c) as a function of J. For 
unligated heme (Fig. 2. a), we do not find any qualitative 
difference in the optical spectra as J is increased from 
OeV to 1.2eV. However, for both FeP(CO) and FeP(0 2 ), 
we find that increasing J gives a double peaked structure 
at uo « 2eV, which is also observed in experiments (see 
inset of the Fig. 3 of the paper). This double peaked 
structure, absent for J = OeV, suggests also that J is im- 
portant to describe the heme molecule. We note that the 
best fit to the experiments is obtained for J = 0.9eV, and 
hence we can furthermore support that the Hund's cou- 
pling is consequent in heme based on phenomenological 
grounds. 
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FIG. 1: (Colors online) a) Iron occupation rid as a function of the Hund's coupling J for different Coulomb repulsion U for 
FeP-d. The data discussed in the paper are reproduced here for comparison (U = 4eV, bold line). The sharp transition in 
the iron density is present for all repulsions, b) Iron occupation rid for two different temperatures T = 100K (diamonds) and 
T = 293K (circles). 




co [eVJ oo [eV] co [eV] 



FIG. 2: Optical spectra for a) FeP-d, b) FeP(CO) and c) FeP(0 2 ) for different values of J (from top to bottom): J = OeV 
(bold line), J = 0.8eV (dotted line), J = 0.9eV (dashed line), J = 1.2eV (full line). The arrows highlight the presence of a 
single peak structure for J — OeV, which splits in a double peaked structure as J increases. For clarity, the curves were shifted 
on the vertical axis. 



EFFECT OF THE HUND'S COUPLING ON THE 
ORBITAL POLARISATION AND BINDING TO 
LIGANDS 

As shown in the Fig. 3a of the paper, the effect of the 
Hund's coupling is to reduce the difference in the FeP 
binding energies to CO and O2 (AAE). In this section 
we give more details regarding the mechanism underlying 
the reduction of AAE induced by J. 

In Fig. 3 we show the energy variation for FeP (CO) 



and FeP(02) as the Hund's coupling J increases. We find 
that while the energy of FeP (CO) is only weakly affected 
by J, there is a sharp drop of the energy of FeP(02) at 
J « 0.8 eV. 

The latter can be understood in terms of orbital polari- 
sations. In particular, we compare in Table I the electron 
occupation of the iron d orbitals for J = and J = 0.8 
eV. We note that the orbital occupations are not integer 
numbers. Indeed, in the presence of hybridisation (see 
Fig. 4.a,b), the atomic iron densities (local observables) 



LU -2 





••FeP(0 2 ) 
•#FeP(CO) 







" 4 o!5 

J[eV] 



FIG. 3: Energy: variation of the total energy AE = 
E(J) - E(J = 0) of FeP(0 2 ) (red circles) and FeP(CO) (blue 
circles) as a function of the Hund's coupling. While there is 
a drop in the energy for FeP(02) at J ~ 0.8, the energy of 
FeP(CO) is only weakly affected by the Hund's coupling. 





J 


d x 2_ y 2 


d 3z 2_ r 2 


d xz 


dec y 




FeP 





0.85 


1.86 


1.24 


1.98 


0.82 


FeP 


0.8 


1.10 


1.75 


1.08 


1.14 


1.08 


FeP(CO) 





1.06 


0.86 


1.99 


1.06 


1.99 


FeP(CO) 


0.8 


1.14 


1.33 


1.16 


1.05 


1.85 


FeP(0 2 ) 





0.72 


1.82 


1.25 


1.87 


1.28 


FeP(0 2 ) 


0.8 


1.03 


1.07 


1.18 


1.97 


1.09 



TABLE I: Average occupations of the iron d orbitals for 
FeP, FeP(CO) and FeP(0 2 ), for J=0 and J=0.8. 



are not conserved quantum numbers, and can take frac- 
tional values, which is a signature of valence fluctuations. 

However, due to the form of the hybridization of the 
iron atom in FeP, the orbitals (d 3z 2_ r 2 ,d xy ) are almost 
filled and hence do not fluctuate (Table I) in the absence 
of the Hund's coupling. We also find that for FeP (CO) 
and FeP(02), the orbital which are almost full are respec- 
tively the (d xz ,d yz ) and the (d 3z 2_ r 2 ,d xy ) orbitals. 

The main effect of the Hund's coupling J is to increase 
the magnetic moment of the iron atom, which in turn 
yields a concomitant reduction of the iron density n^. 
For FeP, there is a reduction of the iron density of 0.52 e 
(see Table III), induced by the build up of the self en- 
ergy at the iron site due to the Hund's coupling and the 
Coulomb repulsion, and a reduction of 0.25 e of the neigh- 
bour N atoms. The charge is transferred to the two hy- 
droxyl chains which are electronegative (each hydroxyl 
chain contains two O atoms). 

We find that in FeP (CO), the orbital mostly affected by 
an increase of J is the d xz : the occupation of the latter 
is reduced from 2 e ( J — 0) down to 1.16 e (J = 0.8eV). 
For FeP(02), the orbital which is strongly affected by J 
is the d 3z 2_ r 2 orbital: the occupation of the latter is re- 
duced from 1.8 e down to 1.07 e. Moreover, we find that 
the reduction of the occupation of d xz in FeP (CO) and 
the reduction of the occupation of d 3z 2_ r 2 in FeP(02) 
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J 


d x 2_ y 2 


d 3z 2_ r 2 


d xz 


dxy 


d yz 


FeP 





0.30 


0.15 


0.28 


0.05 


0.28 


FeP 


0.8 


0.31 


0.20 


0.39 


0.38* 


0.39 


FeP(CO) 





0.37 


0.37 


0.03 


0.48 


0.03 


FeP(CO) 


0.8 


0.37 


0.40 


0.46* 


0.48 


0.19 


FeP(0 2 ) 





0.35 


0.20 


0.34 


0.17 


0.34 


FeP(0 2 ) 


0.8 


0.37 


0.43* 


0.45 


0.08 


0.46 



TABLE II: Internal magnetic moment of the iron d orbitals for 
FeP, FeP(CO) and FeP(0 2 ), for J=0 and J=0.8. The inter- 
nal orbital magnetic moment is obtained by S a = \J (S Q S a ), 
where a is an index for the d orbital and S a is the spin oper- 
ator of the orbital a. The star highlights the orbital with the 
largest moment increase. 



atom 


An{r) 


Iron d orbitals 
Nitrogen ring 
hydroxyl groups 


-0.52 
-0.25 
+0.77 



TABLE III: Variation of the charge An(r) = n(r, J = 0.8) - 
n(r, J = 0) in FeP induced by the Hund's coupling. 

are consistent with a build up of the magnetic moment 
in the latter orbitals (see Table II). We note that both the 
Coulomb repulsion U and the Hund's coupling J have an 
importance here, and promote different many body con- 
figurations. The Coulomb repulsion U tend to suppress 
doubly occupied many body configurations, whereas the 
Hund's coupling tends to generate many body configu- 
rations with aligned spins. In particular, a physical con- 
straint on the Hund's coupling J is J < U (see Ref. [30]). 

The reduction of the charge of the d xz orbital in 
FeP(CO) and of the d 3z 2_ r 2 orbital in FeP(02) is ex- 
pected to reduce the Coulomb energy, which penalises 
doubly occupied orbitals. Moreover, the charge reduc- 
tion is also expected to yield an optimisation of the ki- 
netic energy, since transfer of electrons to a filled or- 
bital are impossible. The latter is however not effec- 
tive in FeP(CO), since the d xz orbital does not hybridise 
strongly (see Fig. 4. a) and hence no drastic change is ex- 
pected in the kinetic energy, but the latter is expected to 
be important for FeP(02) since the d 3z 2_ r 2 orbital hy- 
bridise significantly (see Fig. 4.b). This accounts for the 
further energy reduction in FeP(02) at J w 0.8eV, not 
present in FeP (CO), as observed in Fig. 3. 

QUANTUM DYNAMICS: TRANSIENT 
MAGNETIC RESPONSE SIMULATIONS 

We analyzed the quantum dynamics of heme by 
means of out-of-equilibrium calculations performed us- 
ing a recently-developed implementation of the Keldysh 
formalism adapted for Lanczos DMFT solvers [55]. Our 
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FIG. 4: Iron dynamical hybridization: Imaginary part of 
the hybridization A" as a function of the matsubara frequen- 
cies iujn of the iron d orbitals for a) FeP(CO) and b) FeP(02). 
The e g orbitals have the largest hybridization. Note that the 
hybridization is independent of the many body parameters U 
and J. 
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FIG. 5: Transient magnetic response: Time depen- 
dence of the subspace magnetization wid — (n^ — n^)H/2, after 
an initial quench in a magnetic dipole field, for heme species 
a) FeP-p, b) FeP-d, c) FeP(CO) and d) FeP(0 2 ). 



calculations allow us to directly tackle the issue of dy- 
namical relaxation after a quantum quench; the results 
are shown in Fig. 5. Specifically, we applied a quench 
in molecular spin polarization at t = 0, by perturbing 
the ground state with a magnetic dipole impulse applied 
to the iron 3d impurity subspace and thereafter study- 
ing the transient behavior of the spin polarization. The 
heme molecule plays the role of a dissipation bath, such 
that after a long enough time the system relaxes to the 
equilibrium. The relaxation times for unligated deoxy- 
heme, at J = 0.8 eV, were found to be approximately 
t = 80 fs (Fig. 5.a) and t = 40 fs (Fig. 5.b) for FeP-p 
(unligated, planar) and FeP-d (unligated, domed), re- 
spectively. We find that in the FeP-d system, only the 
d x 2_ y 2 orbital is strongly spin polarized by the quench, 
whereas, in FeP-p, the d 3z 2_ r 2 orbital is also partially 
occupied in the ground-state and hence also participates 
in the magnetic response, yielding transient oscillations 
that are longer-lived overall. Our result suggests that the 
dome-shaped configuration protects heme from external 
perturbations to some extent, reducing the hybridization 
between the d x 2_ y 2 and d 3z 2_ r 2 orbitals by further lifting 
their degeneracy. For FeP(CO) (Fig. 5.c), we observe a 
very long relaxation time (beyond our numerical reach); 
magnetically perturbed FeP(CO) becomes trapped in a 
metastable state with a slowly-decaying internal magne- 
tization. Finally, for FeP(02), we find an intermediate 
relaxation time of t = 80 fs (Fig. 5.d), and transient 
response rather different to that of FeP(CO), a differ- 
ence which appears to be primarily due to a magnetic 
equivalence among the t2 g orbitals of FeP(02) which is 
absent in FeP(CO). These results highlight the surpris- 
ingly long time-scale and system dependence associated 
with quantum spin relaxation. Photo-excitation of car- 
boxyheme, FeP(CO), has previously been discussed in 
the context of both femtosecond spectroscopy [56] and 
time-resolved resonance Raman spectra [57], where the 
presence of excited states with a lifetime of 300 fs was at- 
tributed to iron-to-porphyrin charge-transfer excitations. 
We suggest, based on the results of our out-of-equilibrium 
simulations (Fig. 5), that circularly-polarized two-photon 
spectroscopy may be used as a sensitive probe for the ge- 
ometric configuration and ligation state of heme. 
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